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We present numerical relativity simulations of nine-orbit equal-mass binary neutron star covering 
the quasicircular late inspiral and merger. The extracted gravitational waveforms are analyzed 
for convergence and accuracy. Second order convergence is observed up to contact, i.e. about 3-4 
cycles to merger; error estimates can be made up to this point. The uncertainties on the phase and 
the amplitude are dominated by truncation errors and can be minimized to 0.13 rad and < 1 %, 
respectively, by using several simulations and extrapolating in resolution. In the latter case finite- 
radius extraction uncertainties become a source of error of the same order and have to be taken into 
account. The waveforms are tested against accuracy standards for data analysis. The uncertainties 
on the waveforms are such that accuracy standards are generically not met for signal-to-noise ratios 
relevant for detection, except for some best cases using extrapolation from several runs. A detailed 
analysis of the errors is thus imperative for the use of numerical relativity waveforms from binary 
neutron stars in quantitative studies. The waveforms are compared with the post-Newtonian Taylor 
T4 approximants both for point-particle and including the analytically known tidal corrections. 
The T4 approximants accumulate significant phase differences of 2 rad at contact and 4 rad at 
merger, underestimating the influence of finite size effects. Tidal signatures in the waveforms are 
thus important at least during the last six orbits of the merger process. 



I. INTRODUCTION 

An exciting possibility to reveal and to study the na- 
ture of the neutron star (NS) interior is provided by the 
detection of gravitational waves (GWs) from binary neu- 
tron star (BNS) mergers. The GWs emitted by BNS 
during the late inspiral and merger are sensitive to finite 
size effects, and in particular to the tidal interaction be- 
tween the bodies, thus they are quantitatively dependent 
on the star parameters and, in turn, on the equation of 
state (EoS). Ground-based interferometers are sensitive 
to the last 10 orbits of a typical equal-mass binary sys- 
tem of mass ^2.8 M©, which roughly corresponds to the 
frequency range 400 — 1500 Hz. During this phase tidal 
effects are expected to be significant. 

A general relativistic perturbative theory of tidal in- 
teractions has been developed in recent years [irllj- 
These results have been incorporated into the post- 
Newtonian (PN) formalism Jj-Q; thus permitting the 
extension of phasing formulas to tidally interacting bi- 
naries, as well as into the effective-one-body (EOB) 
model 0- 

An exact and quantitative evaluation of the dynamics 
and of the waveforms during the merger process requires, 
however, the solution of the full nonlinear Einstein equa- 
tions. In particular, numerical relativity (NR) simula- 
tions are to date the only tool to tackle the problem (see 
e.g. [||-[l2j for recent works in the field and [l3|-[l5[ for 
reviews). 

Numerical relativity data have been used in combina- 
tion with PN methods in order to assess the detectability 
of tidal effects and the accuracy of the parameter esti- 
mated from GW measurements [l6j|. A more recent ap- 
plication of NR results concern their use for calibrating 
the tidal-EOB model [|, 03 ■ These works highlight the 



importance of using NR waveforms and analytic results 
in order to quantitatively evaluate the impact of tidal 
effects in the waveforms. 

An aspect which deserves a more detailed assessment 
than is currently available in the literature is the accu- 
racy of the NR waveforms. The estimates of error-bars on 
phase and amplitude of BNS waveforms, as well as their 
assessment against accuracy standard for detection fl8l — 
|20T |. is of fundamental importance for any quantitative 
study. In constrast to binary black holes (BBHs) data, 
whose quality for data analysis purposes is well docu- 
mented, see e.g. [2lTj24| . convergence and uncertainties 
in BNS simulations arc so far poorly investigated. To 
our knowledge, the only analysis of truncation errors 
have been performed in [25j . and, more exhaustively, 
in [2tij . Both works found that the waveforms are sec- 
ond order convergent up to merger but they are limited 
to short runs (three orbits) and do not consider accu- 
racy criteria for detection. In Q the same initial data 
as in this work have been evolved, error bars consider- 
ing finite-extraction effects and truncation errors have 
been estimated but without performing convergence tests 
and using only two simulations. In their conclusions 
the authors stressed the need for a detailed error bud- 
get based on convergence measurements. Since the nu- 
merical treatment of the matter (hydrodynamics) makes 
very challenging to obtain accurate waveforms (in com- 
parison with BBHs simulations), a precise and rigorous 
assessment of their quality is urgent. 

In this work we report results about the accuracy of the 
waveforms extracted from BNS simulations. Focusing on 
an example configuration, we consider nine-orbit simu- 
lations employing different resolutions, and reaching the 
highest resolutions for production runs used so far. We 
discuss convergence in detail, compute error-bars of the 



2 



waveform amplitude and phase, and test them against ac- 
curacy standards for detection. The NR waveforms are 
then contrasted with the post-Newtonian (PN) T4 phas- 
ing formula, both for point-particle and including the 
leading-order and next-to-leading-order tidal corrections 
analytically known. 

The structure of the paper is as follows. In Sec. |TT] 
we review our methodology. In Sec. IIIII the simulated 
binary dynamics is presented. In Sec. IIVI we analyze 
the waveforms. In particular, convergence is analyzed in 
Sec. lIV Al the influence of finite-radius extraction is ana- 
lyzed in Sec. IIVB1 and the waveforms are tested against 
accuracy standards in Sec. lIVD] In Sec.|V]the numerical 
waves are compared with the analytic post-Newtonian 
T4 formula including tidal effects. 

Dimcnsionlcss units c=G=M@=l are used in this pa- 
per, unless otherwise stated. 



II. THEORETICAL AND NUMERICAL SETUP 

In this section we outline the theoretical and numerical 
frameworks employed in this paper, and we describe the 
setup of the simulations. More details on the methodol- 
ogy are given in [26|, [27j and references therein. 

The present study relies on evolutions of BNS initial 
data within 3+1 numerical relativity, using the BSS- 
NOK formulation [2SM30| of Einstein equations coupled 
with the general relativistic hydrodynamics (GRHD) sys- 
tem [3l|. The gauge is specified by the 1+log lapse 
and Gamma-driver-shift [32H35| . using the same expres- 
sions and parameters of Sec. II of [26|]. In particular 
the damping parameter for the shift equation is set to 
r) = 0.3. Gravitational waves are extracted from the 
numerically generated spacetime by using the Newman- 
Penrose scalar, i/> 4 . The projections onto spin weighted 
spherical harmonics, i.e. the multipoles of the radiation, 
are evaluated on extraction sphere of finite coordinate 
radius r. The same conventions of [13] are employed 
here. The diagnostic quantities discussed in the following 
are the ADM mass, Madm, the Hamiltonian constraint, 
Ham, and the rest-mass integral, Mo, Cf. Eq. (31) of I26l| . 

The code employed in this work is the BAM code |26l . 
l27l . l36l [I?]], which implements finite differencing meth- 
ods on Cartesian refined meshes. The evolution algo- 
rithm is based on the method of lines and explicit Runge- 
Kutta methods (third order in this work). A combi- 
nation of centered and lop-sided standard finite differ- 
ences in space is used for the metric fields, see [27|]- In 
this work fourth order operators are employed, together 
with sixth order artificial dissipation. The algorithm 
implemented for the matter is a robust high-rcsolution- 
shock-capturing scheme based on a central scheme for 
the numerical fluxes (3^-(42j. Both the time stepping 
and the spatial refined mesh are shared with the met- 
ric system. The interface fluxes are computed by the 
local Lax-Friedrichs (LLF) central scheme [39|, HO] , while 
reconstruction is performed with the third order convex- 



essentially-non-oscillatory (CENO) interpolation j43l |44| . 
Mesh refinement is provided by a hierarchy of cell- 
centered nested Cartesian grids and Berger-Oliger time 
stepping. Metric variables are interpolated in space by 
means of fourth order Lagrangian polynomials and mat- 
ter conservatives by a fourth order weighted-essentially- 
non-oscillatory (WENO) scheme [45[. Interpolation in 
Berger-Oliger time stepping is performed at second order. 
Some of the mesh refinement levels can be dynamically 
moved and adapted during the time evolution according 
to the technique of "moving boxes", e.g. [27| . 

Initial data are chosen from quasi-cquilibrium con- 
figurations of irrotational equal-mass binaries in quasi- 
circular orbits [ill H3- The configuration selected for 
this work is a binary with ADM mass Madm = M = 
3.00506, rest-mass M = 3.250, and angular momentum 
Jadm = 9.716. The initial proper relative separation is 
d ~ 50 (~ 70 km) corresponding to a GWs frequency 
fo = 0.0019 (394 Hz). The compactness of each star 
in isolation is 0.14. The EoS for the fluid is the poly- 
tropic one, with adiabatic index r = 2. The initial con- 
figuration is computed with a multidomain spectral code 
which solves the Einstein constraint equations under the 
assumption of a conformally flat metric. The code is 
based on the Lorene library [48| and provided by the 
NR group in LUTH (Meudon). These initial data rep- 
resent to date the most accurate computation of equilib- 
rium BNSs and they are publicly available on the web. 
The same initial data were used for the evolutions dis- 
cussed in [t| [I?} ■ 

Evolutions were performed for the r — 2 polytropic 
EoS, i.e. we consider the fluid isentropic and neglect ther- 
mal effects. In (26^ we have shown that, in agreement 
with the physical expectation, waveforms computed with 
both the polytropic EoS and the ideal gas EoS (which in- 
cludes in a rough way thermal effects) are indistinguish- 
able within the simulation errors, at least up to contact, 
while significant differences accumulate during merger 
and the HMNS phase. In this work we are mainly inter- 
ested in the inspiral phase, hence thermal effects can be 
neglected; we will consider thermal effects in the present 
setup in future work. 

Gravitational waves were extracted at levels I = 1, 2, 3. 
Several resolutions and a single grid setup were employed. 
The latter is composed of a fundamental grid level, I = 0, 
and seven refinement levels from I = 1 to Z max = 7; 
four refinement levels are moving, I = 4,5,6,7. The 
only symmetry assumed is reflection symmetry about the 
3 = plane, i.e. the numerical domain is restricted to 
z > 0. The grid configurations, as well as the perfor- 
mances of the runs are reported in Tab. |TJ The grid set- 
tings are similar to those of other codes, e.g. j49[. The 
highest resolution employed here for run HH6 is slightly 
(3%) higher than the maximum resolution used to date 
on BNS simulations employing mesh-refined-Cartesian- 
grid-based codes @. All the runs were performed with 
Courant-Friedrich-Lewy (CFL) factor of 0.25. The (self) 
convergent series are formed by triplets of runs, that 
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TABLE I: Summary of the grid configurations and of the runs. Columns: name of the configuration, maximum refinement 
level, minimum moving level, number of points per direction in the moving levels, resolution per direction in the level / = Z max , 
number of points per direction in the nonmoving levels, resolution per direction in the level / = 0, number of processors, maximal 
memory usage, and average speed in term of the mass of the configuration evolved (M = 2.998 Mq) including checkpointing 
and initialization (reference machine: JUROPA cluster). 
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in the following will be denoted as HH{LMH}, where 
L, M, H correspond to the low, medium, and high res- 
olution employed. Self-convergence tests can be biased 
by the choice of the resolutions employed. An "opti- 
mal" setup would require that (i) the ratios between the 
low and medium and medium and high resolutions are 
hh/hM — hM/hjj — 2; and (ii) the scaling factor is at 
least of order two, SF = {h r L - h r M )/(h r M - h r H ) > 2, 
where r is the convergence rate. It is difficult to obtain 
an optimal convergent series since low resolutions are too 
inaccurate and differ even in a qualitative way [26| . Con- 
sidering the criteria above, the "best" convergent series 
is HH{235} or HH{236}. 



III. OVERVIEW OF THE BINARY DYNAMICS 

In this section we summarize the binary dynamics and 
present some diagnostic of the simulations. 

The binary evolves for about nine orbits dynamics be- 
fore merger, when a hyper-massive-neutron-star (HMNS) 
is formed. The latter oscillates non linearly in time, 
loses angular momentum by GW emission increasing its 
compactness, and, finally, collapses to a black hole sur- 
rounded by a disk rapidly accreting. 

The evolution of the maximum rest-mass density is re- 
ported in Fig. Q] (left), together with the proper distance 
(right). During the inspiral the maximum of the rest- 
mass remains constant as expected; the proper distance 
shows some residual eccentricity from the initial data. 
The latter is larger during the first three orbits and then 
progressively radiated away, although not completely. 
The stars touch each other about 1.5 orbits before merger 
(t/M > 2050), which happen at t m /M = 2259 for run 
HH6 (see Sec. IIVI for the definition of merger used also 
in this work). After the merger the maximum of the 
rest-mass density increases indicating the compactness 
of the HMNS increases; it reaches a peak during the 
collapse then drops down to the densities of the accre- 
tion disk. The quasi-radial oscillations of the HMNS are 
also visible before the collapse. An apparent horizon is 
formed at tp^/M ~ 2475 (run HH6), the mass and spin 
of the final puncture describing the black hole [H^-HH are 
M BH = 2.955 ± 0.005 and a BH = 0.80 ± 0.01. The latter 



values are computed from the irreducible mass and spin, 
after an initial transient in the BH formation. 

Overall the new simulations at higher resolution con- 
firm our previous findings about the merger outcome (2fj| | : 
the HMNS experiences a delayed collapse (8(| while a 
prompt collapse seems an artifact of lower resolution runs 
(see HH2). Note also that the use of lower resolutions re- 
sults in earlier mergers. 

The evolution of the rest-mass is reported in 
Fig. [21 During the inspiral it is conserved up to 
max(AMo/Mo) < 1% for runs HH4 and higher resolu- 
tions. At the collapse it drops several order of magni- 
tude, similarly to the maximum of rest-mass density. As 
explained in [5(J this effect is produced by the gauge con- 
ditions which, handling the singularity formation, stretch 
the numerical grid effectively moving grid points to larger 
proper radii. The values of Mq at late times are an es- 
timate (upper limit) for the rest-mass of the accretion 
disk, which is below 1% of the initial rest-mass. Note the 
strong dependence of the result on the resolution. 

The convergence of the L2 norm of the Hamiltonian 
constraint is reported in Fig. [3] The different data set 
are rescaled for second order convergence to the highest 
resolution one. As one can observe from the figure, the 
lines are superposed during the inspiral while they pro- 
gressively differ from the contact, towards the merger. 
After the merger convergence is not measurable. 

Finally we comment about the ADM mass conser- 
vation. We computed finite-radius approximations to 
the ADM mass, MadmW, by integrals over coordinate 
spheres as in the case of the GW, considering two differ- 
ent formulas: (a) Eq. (54) of 27| . and (b) the integral of 
the conformal factor only. The ADM mass is defined for 
the limit of large spheres, r — > oo. The value at finite r is 
a coordinate dependent quantity, and, for large but finite 
r, it suffers of resolution problems in the outer levels. In 
our setup the calculation is not accurate enough to make 
quantitative statements and extrapolation r — > oo does 
not seem to improve the results. We observe anyway 
a consistency between Madm(0 and the energy of the 
emitted GW within the 1 % level. 
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FIG. 1: BNS dynamics, (left) Evolution of the maximum of the rest-mass density for different resolutions, (right) Evolution 
of the proper separation between the two stars. 
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FIG. 2: Evolution of the rest-mass density for different res- 
olutions. 
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FIG. 3: Convergence of the L2 norm of the Hamiltonian con- 
straint. The data set are rescaled for second order convergence 
to the highest resolution one. 



IV. WAVEFORMS 



The total gravitational energy radiated during the 
merger process is about 1 % of the initial ADM mass. 
About 99 % of the energy radiated during the inspiral is 
emitted into the (£, m) = (2, 2) multipolar channel: the 
latter is also responsible for about the 92 % of the energy 
emitted during the whole simulation. In the following we 
will consider only the (2,2) mode. Figure 0] (left) shows 
the real part, the imaginary part, and the absolute value 
of the GW multipole r h^i extracted at r = 750 (250 M) 
and from the HH6 run. All the plots relative to the wave- 
forms are in term of the retarded time without changing 



notation. We use thus t — > t — r*, where the tortoise 
radius is computed as r» = R + 2Mlog(i?/(2M) - 1), 
with R(r) the Schwarzschild radius corresponding to the 
coordinate extraction (isotropic) radius r. 

The waveform is characterized by the chirp-like shape 
typical of the quasi-circular inspirals, after about 18 cy- 
cles it peaks, and then shows a more complicated struc- 
ture with multiple maxima in amplitude and progres- 
sively higher frequencies. We formally define the merger 
time, t m , as the time corresponding to the peak of the 
amplitude of r hyi 0, H(| . The signal after the merger 
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is characterized by the emission from the HMNS. We re- 
ported an analysis in [26|, see also [53| for recent work. 
More details will be given in a future work, in the follow- 
ing we will focus on the inspiral waveforms. 

The simulation output is the multipolc of the curva- 
ture scalar, tfj 22 . The actual GW strain, h, is recovered 
from ^ 4 by integrating the relation h — tp 4 . The integra- 
tion is not straightforward in the case of nois y n umerical 
data. We employ the method described in [54j, i.e. we 
perform the integration in the Fourier frequency domain 
by applying a fixed-frequency-high-pass filter, following 
closely [261 ], The cutting frequency used in this work is 
/cut = 0.0016 < Jo- As shown in Fig. [U the waveform is 
affected by some amplitude modulations mostly present 
at early times. Their origin may be due to the residual 
eccentricity contained in the initial data. 

In the following waveforms will be split into phase and 
amplitude according to the notation, 

r h 22 = A 22 exp (-i$) , r ip 22 = °22 exp (— i<j)) ■ (1) 

The instantaneous GW frequency is lo — — 3(/i//i), and 
is plotted in Fig. [4] (right) in case of the (2, 2) multi- 
pole. It increases monotonically during the inspiral, and 
reaches the value M u>22(t m ) ~ 0.123 at the merger. At 
t/M ~ 2400 it shows a signature of the HMNS, and at 
later times increases to the quasi normal modes (QNMs) 
frequencies of the final black hole (BH). Note that the 
GW frequency drops to zero at t/M ~ 2300, corre- 
sponding to a minimum of the amplitude and to a quasi 
spherical shape of the stars (26|. The frequency of the 
BH fundamental QNMs can be extracted from the GW 
frequency, however a cleaner equivalent signal is pro- 
vided by the frequency of r ip 22 ■ We found for run HH6 
/qnm — 6.47 kHz {Mw 22 ~ 0.6), in 2 % agreement with 
the estimate obtained from the horizon quantities. 

In the following the accuracy of the numerical wave- 
forms is assessed. We stress that here, for the first time, 
phase and amplitude errors are measured precisely and 
consistently from convergence tests. 

A. Convergence 

In this section we present the results concerning the 
self-convergence of the inspiral waveforms. The conver- 
gence series HH{236} is discussed as an example, similar 
results are obtained for HH{235}. We focus on the ex- 
traction radius r = 750 and on rijj 22 . Similar results are 
found for r h 22 . 

In Fig. [5] the self-convergence test is shown. The scal- 
ing factor is SF(2) = 1.8. The differences are noisy so the 
figure employs a standard Savitzky-Golay averaging filter 
for a better visualization; results are not affected anyway. 
The waveforms show compatibility with second order 
self-convergence during the inspiral up to t/M ~ 2000. 
At later times they become, as other quantities, progres- 
sively over-convergent, and after the merger the conver- 
gence order can not be established. We observe here 



a common finding in NR simulations (e.g. (2f| HH, HBj]): 
as long as only the bulk motion of the matter is im- 
portant, the numerical methods employed do quite well 
modelling the inspiral due to GW emission, but degrade 
when strong field and matter dynamics develop. 

The over-convergence behavior appearing at late times 
is probably due to the run HH2, but also to the fact 
that, when the stars come in contact, the effective order 
of the (nonlinear) numerical scheme for hydrodynamics 
probably drops below the second order (in norm). As in 
previous shorter runs (26[, the phase is not exactly con- 
vergent at rate two but at lower rate (between one and 
two). Note that, differently from previous works on BNS 
and consistently with (26j , we do not align the waveforms 
for the convergence tests. The gravitational energy car- 
ried by the (2, 2) mode also shows approximately second 
order self-convergence. 

The interpretation of these data can be delicate be- 
cause several sources of systematic errors are not com- 
pletely under control: the exact expected convergence 
rate, the role of different grid setup, the limited and 
not optimal choices of resolutions for the convergent se- 
ries, etc. Our findings, however, appear consistent and 
sufficiently robust; the second order rate is expected, in 
convergence regime, by basic arguments, and the diag- 
nostic quantities of Sec. IIIII show second order conver- 
gence in norm. The results seem to indicate that second 
order convergence can be confidently assumed up to con- 
tact, or, equivalently, to M oj 22 = 0.07. In the following 
we will assume second order convergence for the extrapo- 
lation of the inspiral waveforms up to merger, errors will 
be given both for M u>22 _ 0.07 and for M u>22 _ 0.1. 
The reliability of the latter estimate is not clear. 



B. finite-radius extraction 

In this section we study the uncertainties on phase 
and amplitude related to the computation of wave- 
forms at finite-extraction radii. We consider several 
extraction radii r = 200, 300, 400, 500, 750 (or R ~ 
203, 303, 403, 503, 753) from run HH6 and r^ 2 . 

The differences in amplitude and phase extracted at 
a given radius with the previous, e.g. A*(f> 22 (Ri) = 
4>22(Ri)~4>22(Ri-i), are shown in Fig. [5] Both amplitude 
and phase increase for higher extraction radii. The dif- 
ferences are bigger at earlier times, the phase differences 
at early times scale approximately as 1/r, while ampli- 
tude differences approximately as 1/r 2 . For radii r > 400 
differences in amplitude are < 2% and they seem to satu- 
rate. By contrast differences in phase keep on increasing 
and between r — 750 and r = 500 they are ~ 0.1 rad. 
The differences become progressively smaller towards the 
merger. 

Following previous works (56l46ll | , an approximation of 
the waves at null-infinity can be obtained by simple 1/R- 
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FIG. 4: r/122 waveform, (left) Real part, imaginary part, and amplitude; (right) instantaneous GW frequency. Extraction 
radius r = 750, run HH6. 




extrapolation, 

K 

F(t,R)=Y. F »^ R ~ k > ( 2 ) 

fc=0 

where F(t, R) is either the phase or the amplitude, t is 
the retarded time, and F (t) is the extrapolated value. 
In I(i2l (i 1 the robustness of the extrapolation proce- 
dure has been assessed against null-infinity waveforms 
from equal-masses BBH inspirals computed with the 
Cauch y-ch aracteristic extraction (CCE) method [65l - l67j . 
In [68, 69], by mean of BH perturbation theory on hyper- 
boloidal foliations, it has been shown that, in case of an 
unambiguous definition for the background and for the 



retarded time, the extrapolation reproduce null-infinity 
waveforms up to their numerical uncertainties for enough 
high values of K > 3. 

Figure [7] shows the differences between the extrap- 
olated value for different K and the reference radius 
r = 750. Note that waveforms at different radii are 
not shifted in time or phase, but only considered against 
the retarded time. The fit errors, SF, are computed 
at the 68 % confidence level, and distributed quite uni- 
formly in the inspiral. Hence, their average, (SF), can 
be use as a meaningful measure of the fit quality in com- 
parison with the differences, AF, between the extrap- 
olated waves and the finite-radius extracted ones. In 
case of a linear extrapolation, K = 1, the average fit 
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FIG. 6: Differences between amplitude (left) and phase (right) of waves extracted at successive radii r = 200, 300, 400, 500, 750. 
Run HH6. 



errors are (^22) ~ 0.02 rad and (60,22/0.22) ~ 3 %, 
and the maximum differences with the last radius are 
Aa22/a22 ~ +5 % and maxA022 ~ +0.23 rad. For 
K = 2, the average fit errors are (6cp 22 ) ~ 0.04 rad and 
(5022/022) ~ 3 %, and maximum differences with the 
last radius are maxAa22/fl22 ~ +1% and maxA(/>22 ~ 
+0.21 rad. The use of K > 2 results in more noisy 
data as shown by the figure, and also the fit errors in- 
crease. The "best" extrapolation is thus given by K = 1 
or K = 2. Note however that the fit averaged uncertain- 
ties are about 10 % of the phase difference with the last 
resolution also in the best cases, and that, within this 
uncertainty, both the extrapolations basically agree (see 
right-hand panel of Fig. E}. 

A similar behavior has been observed for the extrapo- 
lation of r /122 ■ In this case however data are less noisy 
and the fit errors are smaller. Specifically we found 
(5^22) ~ 0.006 rad and (5A 22 /A 22 ) < 2 % for K = 1 and 
(<J$) 22 ~ 0.002 rad and (SA22/A22) < 1 % for K = 2. 
The latter is thus preferable. The differences with the 
last extraction radius are reported in Fig. [5J and compa- 
rable in absolute size to those of Fig. [7] 



C. Truncation errors 

In this section we quantify the truncation errors in 
the inspiral waveforms. Richardson extrapolation is em- 
ployed using different data sets and assuming second or- 
der convergence. Extrapolation series are indicated with 
the same notation of convergence series, e.g. HH{23456}. 
Errors are computed as differences with the highest res- 
olution data. We stress that this is a common but opti- 
mistic choice. Also, to avoid underestimates of the errors, 
the whole convergent series is, at least, used in the ex- 
trapolation. 



Figure O shows the differences in amplitude and phase 
between the extrapolated data from HH{23456} and 
those from run HH6. Similar plots were produced for 
r ^22 an d f° r different extrapolation series. The differ- 
ences are negative and increase towards the merger. Be- 
fore the merger the trend changes and they rapidly in- 
crease to positive values. From the argument given at the 
end of Sec. II V A[ we do not expect to have a fully reliable 
extrapolation at the merger, thus proper error estimates 
must be restricted to slightly before that point. In case of 
other extrapolation series the errors are bigger but qual- 
itatively they show the same behavior. The maximum 
absolute errors observed before the merger are reported 
in Tab. |TT] for different extrapolation scries. They are 
computed within the GW frequency intervals M oj 22 = 
[0.0358,0.07] (i.e. / e [/o,/ m «] = [0.0019,0.0037]), 
where the extrapolation is reliable, and also within the 
GW frequency intervals M u; 22 = [0.0358, 0.1]. In the lat- 
ter case they roughly correspond to the minima in Fig. 

As shown by the table, it is necessary to include at least 
four resolutions to obtain a phase error A$22 J$ 1 fad and 
an amplitude error of A^422/^22 < 1 % for M oj 22 < 0.07. 
In this case truncation errors become of the same order 
of magnitude of the finite-extraction effect. The error es- 
timate up to M UJ22 = 0.1 indicate how dramatically the 
errors increase up to merger. This analysis suggests that 
truncation errors represent the main source of uncertain- 
ties in BNS simulations. 



D. Accuracy 

In this section we test the inspiral waveforms against 
accuracy standards for data analysis. As a measure of 
the accuracy we employ the square of the inaccuracy 
functional, I 2 , whose definition is discussed in detail 
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TABLE II: Maximum differences between extrapolated values in resolution and the highest resolution data (run HH6) in 
different quantities during the inspiral. The first three rows refer to the maximum errors for M U22 < 0.07, while the last three 
rows refer to the maximum errors for M LO22 < 0.1. 



Runs 


Aa 22 /a 22 


[%] |A0 22 | [rod] 


|Ayl 22 /A 22 | [%} 


|A$ 22 | [ 


rad] |Acj 22 /oj 22 [%] 


HH{236} 


38 


1.8 


8 


2 


9 


HH{2346} 


5 


0.4 


1 


0.4 


2 


HH{23456} 


1 


0.13 


0.2 


0.13 


0.6 


HH{236} 


MOO 


7 


20 


5 


>100 


HH{2346} 


70 


1.4 
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1.4 


15 


HH{23456} 


13 
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FIG. 7: Differences between amplitude (left) and phase (right) 
of rip22 extrapolated up to different orders, K, and at the 
reference radius r = 750. Run HH6. 
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in Appendix [2] Eq. (|A3|) . Accuracy requirements are 
set to minimal levels and an ideal detectors is assumed. 
As mentioned in Appendix [21 two waveforms are distin- 
guishable depending on the signal-to-noise ratio (SNR, g) 
of the detection. Given a difference between two wave- 



forms (the error bars, in our case), Sh, there always exists 
a sufficiently high SNR such that the difference is signif- 
icant (in our case, the waveforms are inaccurate). The 
point is thus to assess the accuracy with respect SNR 
that are high enough but also realistic for the future de- 
tections. We recall that for equal-masses BBH waveforms 
accuracy standards are achieved for relevant SNR, and 
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FIG. 9: Differences between extrapolated values from series 
HH{23456} and the data at various resolution. Amplitude 
(left) and phase (right) of r Y121 (bottom). Extraction radius 
r = 750. 



waveforms can be considered faithful, e.g. [2ll.l23l.l24j.l70l]. 

By contrast, such analysis for BNS has never been con- 
sidered before. 

We recall that the inaccuracy functional has been 
rarely employed in NR literature, see e.g. (24[, but it 
provides an equivalent measure to the most common mis- 
match functional, M. (see again Appendix fX] for the def- 
inition). The inaccuracy functional is here preferred be- 
cause its value does not depend on the distance between 
the detector and the source or on the normalization of 
the PSD of the detector noise. In the accuracy standards 
the dependency on the distance from the source is then 
moved to the right-hand-side of the expressions. All the 
results presented here can be translated in term of the 
latter considering that t 1 j q m 2M [23l. FfH] . 

In Tab. IIIII we report, for several detector configura- 
tions, the SNR, computed assuming the source at an ef- 



fective distance of 100 Mpc, and the inaccuracy func- 
tional, computed for different waveforms and choices 
of the error bars. Specifically, l 2 [h m ,h x ] is com- 
puted employing as exact waveform (/i x ) three wave- 
form extrapolated in resolution, HH{236}, HH{2346}, 
and HH{ 23456}, and one extrapolated in resolution and 
radius using the series HH{23456} and the K = 2 series 
of Sec. IIVBI The choice of the model waveform, h m , in 
I 2 basically determine the size of the uncertainties. In 
the table we employ the waveform from the highest res- 
olution run HH6 and extracted at r — 750, except for 
the last column which employed the extrapolated in ra- 
dius form the same run. Other choices were considered 
but they are not shown in the table. The Wiener scalar 
product is computed on the interval / <G [/o, /max] a s in 
Sec. IIV Cl Extrapolation in radius is marked with *. 

As already clear from the analysis of Sec. IIV CI the 
inaccuracy decreases by increasing the number of runs 
used in the extrapolation in resolution. The inaccuracy 
further increases if finite-extraction effect are included. 
Let us consider the criteria in Eq. (| A6|) , and the minimal 
requirements e = 0.5 0, [ill] and £m = 0.005 or £m = 
0.035 [lijj . Note that they correspond to mismatches of 
0.5 % and 3.5 %, respectively, where a 3.5 % mismatch 
indicates that no more then 10 % of the signals are lost. 
Overall our results indicate that: 

(i) the extrapolated waveform form the scries 
HH{ 23456} is effectual and faithful for SNR g < 10 for 
most of the configurations if errors are computed from 
run HH6 and finite extraction effects are neglected or 
included in both h x and /i m at the same time; 

(ii) the extrapolated waveform HH{2346} and the ex- 
trapolated HH{23456}* with errors from HH6, are effec- 
tual only for the less restrictive requirement, £m = 0.035, 
and faithful for SNR g< 3; 

(iii) if waveforms arc extrapolated from less runs 
and/or computing errors from runs at lower resolutions 
then HH5, the inaccuracy had always larger or compara- 
ble values to those reported in the table for HH{236}. 

In conclusion, minimal requirements for data analy- 
sis are met if waves are extrapolated in resolution from 
more then four runs and certain optimistic choices for 
the error bars are made. In the other cases waveforms 
are inaccurate. Similar statements can be made if the 
inaccuracy functional is computed up to M U22 = 0.1, 
while obviously its values increase slightly. 



V. COMPARISON WITH POST-NEWTONIAN 
T4 PHASING FORMULA 

In this section we perform a comparison between the 
NR waveforms and PN approximants. The main goal is 
to quantify their agreement/disagreement and the rel- 
ative signature of the tidal interactions on the waves 
during the last nine orbits of the merger process. The 
comparison presented here is not exhaustive; a system- 
atic investigation of the different phasing formulas (see 
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TABLE III: Inaccuracy functional for several configurations. Columns: detector configuration with reference for the noise 
curve, SNR at 100 Mpc (assuming the detector and the binary are optimally aligned), inaccuracy functional for different choice 
of the waveforms. HH{2 — 6} indicates the waveform has been extrapolated in resolution with those runs, * indicates the 
extrapolation in radius is performed. The model waveform, h m , in the inaccuracy functional is always the one from run HH6 
at r — 750, except for the last column where the extrapolated in radius from the same run is used. 





Q 




I 2 


h m , h x ] 










HH6, HH{236} 


HH6, HH{2346} 


HH6, HH{23456} 


HH6, HH{23456}* 


HH6*, HH{23456}* 


advLIGO [72] 


3.6 


0.333 


0.117 


0.043 


0.149 


0.045 


advLIGO NSNS Opt [73] 


5.4 


0.352 


0.126 


0.046 


0.144 


0.048 


advLIGO Narrow Band [73] 


3.2 


0.462 


0.198 


0.072 


0.112 


0.074 


advLIGO High Sens. [73] 


5.1 


0.401 


0.153 


0.055 


0.132 


0.058 


advVIRGO [72] 


5.1 


0.445 


0.182 


0.066 


0.117 


0.068 


ET [74,75] 


52.0 


0.383 


0.142 


0.051 


0.138 


0.054 



e.g. [56|) and fitting models, as well as the investigation 
of different comparison procedures (e.g. [z3), is beyond 
the scope of this work. 

Here we will focus only on the so-called T4 for- 
mula [H, l78l-[82l]. T4pp hereafter, accurate at 3.5 PN 
level. In addition to the point-particle T4, we will con- 
sider a "tidal" T4, T4td hereafter, as proposed in 0- 
0,0. T4td includes the leading-order (LO) and next-to- 
lcading-order (NLO) tidal PN corrections in the dynam- 
ics and the leading-order corrections in the waveform [87| ■ 

A comparison with the T4 approximants and NR data 
have been already considered in 0, Q3] , to which we also 
refer for the precise equations used in this work. The 
analysis there is performed in frequency domain consider- 
ing a certain measure of the phase acceleration, that has 
the advantage of being independent on time and phase 
shifts (and a simple physical interpretation, sec discus- 
sion in Sec. IV) but the drawback of requiring fits of 
the numerical data and a certain fine tuning. The result 
obtained is that T4td NLO accumulates about 2.25 rad 
on the frequency interval M lj-22 € [0.043, 0.057] for the 
model employed here, and about 2tt rad on the same 
interval for a binary with less compact stars. 

In the following both T4pp and T4td are considered 
in a time-domain comparison with the NR waveform. In 
order to be contrasted, the waveforms must be aligned in 
time and phase. Let us make some general comments on 
this point. There is no unique way to align waveforms 
for such comparison, in the literature several methods are 
proposed, e.g. [lfl [53, UM- A priori none of them is free 
from ambiguities or clearly preferable. The alignment 
region is typically chosen after the "initial transient" (or 
adjustment) of the numerical waveforms. The transient 
is related to the use of conformally flat initial data, and 
the main effect (but in principle not the only one) is 
the well-known burst of radiation at early times of the 
simulation [88}. The transient is quite rapid, typically 
within the first orbit, after that the system relaxes to 
the expected quasicircular state [83j . The allowed align- 
ment region is constrained by the validity of the post- 
Newtonian approximation and the length of NR wave- 
forms. The lowest frequency interval, compatible with 



the NR data available and the comment above, may be 
thus preferable. For a quantitative analysis on the length 
requirement of NR waveforms in the BBHs case see [23[ . 

Guided by these considerations, we chose the following 
strategy: (i) NR and PN waveforms are aligned in phase 
and time by considering an interval, [t\ , £2] , in the first 
half of the numerical signal available, where the frequen- 
cies are closer to those of validity of the PN method; (ii) 
following [58| the time shift, A s t, and the phase shift, 
A S <I>, are determine by minimizing the functional, 

G[A s i, A s $] = I ' dt [$(t) - $ PN (i - A s t) - A s $] 2 . 

(3) 

The PN waveforms are then matched to NR ones by ap- 
plying the shifting; (iii) different results are obtained if 
the center and the length of the alignment interval are 
varied. However we observed that the main dependence 
is on the position of the center rather than in the inter- 
val length. For simplicity, we fixed the interval length as 
100 M and vary the position of the center t c e [0, 900] M; 
(iv) The best value of t c is estimated by minimizing the 
mismatch between the PN and NR waveform. 

As a case study we focus on our best NR waveform ex- 
trapolated only in resolution, i.e. HH{23456}. Neglecting 
finite-extraction uncertainties does not particularly af- 
fect the conclusions, also because they decrease towards 
merger. As discussed extensively in ffl\ the error esti- 
mates of Sec. II VI based on the convergence analysis may 
not be the optimal ones to be used in PN comparison. If 
only a certain range of frequencies of the NR waveform is 
of interest, a phase error estimated on that range (i.e. by 
shifting in some way NR waves from different runs) may 
be less conservative and thus preferable for the specific 
application. However, such error estimates suffer of am- 
biguities related to the alignment procedure and we pre- 
fer not to pursue that method. For the purpose of this 
section is sufficient and justified to use the errors esti- 
mated in Sec. [IV] that can be, eventually, considered as 
an upper bound to the actual errors (see below). 

Figure [10] shows the real part of the aligned waveforms 
(error bars are not shown there for clarity), and the GW 
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frequencies with error bars. The alignment interval used 
for the analysis in the figure is \t\,tq\jM = [450,550], 
which minimizes the mismatch between the NR and T4td 
waveform as described above. Figure fTTI shows the phase 
differences between the waveforms in time (left) and fre- 
quency (right) domain. The PN waveform maintains 
a good phasing for few GW cycles after the alignment 
region and up to t/M ~ 1200. At later times, phase 
differences with respect the PN evolution become posi- 
tive and significant, the largest difference is the one with 
T4pp. At higher GW frequencies than M U22 ~ 0.05 
the PN approximant significantly differ from NR waves, 
and the T4pp rapidly accumulates a phase difference of 
A$22 ~ 2 rad at M w 22 ~ 0.07 and of A$ 22 ~ 4 rad 
at M W22 ~ 0.1. The T4td performs slightly better then 
T4pp at later times, but (not yet calculated) higher or- 
der tidal corrections are important [89j]. This is the cen- 
tral observation here: tidal interactions in the nonlinear 
regime dominate the dynamics and the GW emission at 
least during the last 5-6 orbits of the merger process. A 
similar conclusion can be drawn considering the wave- 
form HH{2346}, but not for HH{236}. In the latter case 
the PN and NR signals are indistiguishable due to larger 
error bars. 

As mentioned above, we did not perform a systematic 
study of different alignment procedures, but a certain 
dependence on the alignment interval was expected and 
observed. The phase differences given above are lower 
bounds, since they are determine by a minimization of 
the mismatch functional. In particular, dropping the step 
(iv) in the procedure outlined above and varying t c <E 
[0, 900] M, we estimated also an upper bound. The latter 
corresponds to an alignment interval centered at t c ~ 
70 M and the accumulated phase are A$22 ~ 3 rad at 
M0J22 ~ 0.07 and of A$ 2 2 ~ 6 rad at M w 22 ~ 0.1. 



VI. CONCLUSIONS 

In this paper we have presented results about the accu- 
racy of NR waveforms from BNS mergers and their com- 
parison with PN methods. The simulations cover nine 
orbits of the late inspiral and the merger phase, they are 
the longest and most accurate BNS simulations to date, 
in terms of the resolution employed and the number of 
runs performed for a single initial configuration. The 
convergence of the waveforms and their uncertainties re- 
lated to truncation errors and finite-radius extraction are 
discussed. For the first time in case of BNS merger wave- 
forms, the accuracy standards for detection have been 
evaluated. The aim of the study is to assess the quality 
of NR waveforms in view of their future use to under- 
stand the physics of the merger and of tidal interactions 
or for data analysis purposes. As a first step in this direc- 
tion, a comparison with the PN T4 waveforms has been 
presented. 

NR waveforms are found to be convergent at second or- 
der rate during the inspiral and up to contact, i.e. until 
the last 1.5 orbits or, cquivalently, for the GW frequen- 
cies Mw 2 2 < 0.07. Later an over- convergent behavior 
is observed, likely due to the numerical treatment of the 
matter and possibly also to the lowest resolution run em- 
ployed in the self-convergence test. The uncertainties on 
the inspiral waves have been estimated by using Richard- 
son extrapolation of different data sets and assuming the 
observed convergence rate where it is valid, finite-radius 
extraction affects were investigated by extrapolating the 
waves to null-infinity. Truncation errors increase towards 
the merger, when the amplitude of the GW become max- 
imum. The maximum errors in phase and amplitude ob- 
served are reported in Tab.HIlfor different extrapolations 
series. The errors related to finite-radius extraction de- 
crease towards the merger and they are generically the 



12 



3 3h 



O 2 

< 

1 

-1 



-NR-T4pp 
NR-T4td 



I 



500 



1000 
t/M 



1500 



2000 



"D 

CO 



CNJ 3- 



< 





- NR-T4pp 






NR-T4td 

















0.04 0.05 



0.06 „ 0.07 

MM 22 



0.08 0.09 



0.1 



FIG. 11: Phase differences T4pp/T4td-NR in time (left) and 
frequency domain (right). The shaded red area at early time is 
the alignment region. The shaded blue area is the uncertainty 
of NR data. 



smaller then truncation errors, but of the same order of 
magnitude of truncation errors in some relevant cases of 
waves extrapolated in resolution. Some of these results 
are compatible with the findings of . 

Accuracy standards have been evaluated using noise 
curves of ground-based detectors, and assuming minimal 
requirements. Results are reported in Tab. IIIII Consid- 
ering the most optimistic error-bars, extrapolated wave- 
forms from five runs are effectual and faithful for de- 
tection with SNR q < 10 for most of the configura- 
tions considered. Considering instead more conservative 
error-bars, or extrapolation in resolution with fewer runs, 
waveforms are neither effectual nor faithful for relevant 
SNRs. 

These facts may affect some of the conclusions of pre- 
vious works where errors of the NR waveforms were (not 



available and thus) neglected, but statements about the 
detectability of (small) effects (EoS, magnetic fields) were 
made. Our results should be taken into account and/or 
reproduced in future works employing NR waveforms for 
data analysis purposes or for comparison with analytic 
methods. 

The NR data have been compared with the prediction 
of the PN T4 formula, both for point-particle (T4pp) 
and including all the analytically known tidal correc- 
tions (T4td). The comparison between the NR and the 
T4 waves has been carried out by aligning them in time 
and phase at low frequencies, and looking at the accu- 
mulated phase difference. The aligned T4pp waveform 
accumulates rapidly a significant dephasing of A 22 ~ 2 
at Mw 22 ~ 0.07, and A0 22 ~ 4 at Mw 22 ~ 0.1. These 
values can be considered lower bounds since in the com- 
parison waves are aligned in such a way to minimize the 
mismatch and error-bars from the convergence test are 
employed. The inclusion of tidal corrections does not 
reduce the phase difference more than a fraction of a ra- 
diant. The results suggest that tidal interactions are very 
amplified in a strong field and nonlinear regime and play 
a significant role already during the last nine orbits. As 
already observed in Q the analytically known LO and 
NLO tidal terms in the T4 PN approximant are not suf- 
ficient to match the NR waveform. 

In summary, our work indicates that NR waveforms 
from BNS are physically "reliable" because convergent 
and comparable with PN at sufficiently low frequencies. 
The measured uncertainties are such that the NR wave- 
forms from BNS may not be sufficiently accurate for data 
analysis purposes, unless data extrapolated from several 
runs are employed. For data analysis applications, an 
error estimate based on relative as well as absolute com- 
parisons (aligning/not-aligning the waveforms) will be 
relevant. A very careful evaluation of the waveform un- 
certainties is unavoidable for their use in quantitative 
studies. 

Future work will be devoted to a more comprehen- 
sive comparison between the NR waves and analytic PN 
fitting models, either in the standard Taylor form or re- 
summed one (EOB), to extend these results to different 
mass ratios and EoSs, and to the investigation of different 
grid setup and higher order methods for the treatment of 
the matter. 



Appendix A: Accuracy standards 

In this appendix the accurac y st andards used in this 
paper are discussed. We follow [US III H ■ 

Given two real time series (waveforms), h^ m (t), and 
their Fourier complex transform, ft. x ,m(/), the Wiener 
scalar product is defined as, 



(h x \h m ) = 45ft / df 



Mf)h* m (f) 

Sn(f) 



(Al) 
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where S n (f) is the one-sided power spectral density of 
the detector noise. The (squared of the) norm associated 
to the Wiener product is ||/i x ||w = (^x|^x)- The Wiener 
product is real and symmetric, the associated norm is 
positive definite. The Wiener product provides a measure 
in the waveform space [§3 ] . The mismatch functional is 
defined as, 



M[h m ,h x ] 



{K\K 



l^x||v 



(A2) 



and it is often employed to discuss accuracy standards 
for BBH NR waveforms, see e.g. plM^ . In this work 
we will mainly focus on (the square of) the inaccuracy 
functional |76| . defined as 



[^m, ^x] 



h x \ 



(A3) 



The inaccuracy functional has been considered for NR 
waveforms, for example, in [24^ . and provides an equiva- 
lent measure to the mismatch functional. 

Assuming an ideal detector (i.e. neglecting calibration 
errors) , a waveform h m (the "model" ) is indistinguishable 
from the waveform /i x (the "exact" ) , if and only if their 
difference dh = h nl — /i x satisfies, 



H^Hw < 1 



(A4) 



Eq. (|A4j) represents an accuracy requirement for mea- 
surement purposes, and it determines the faithfulness of 
the model waveform. A less restrictive requirement can 
be given for detection purposes, and it is related to the 
effectualness of the model waveform. A sufficient condi- 
tion is, 



\\Sh\\ w < ^/2e M Q , 



(A5) 



where g = \\h x \\ w is the optimal signal-to- noise ratio 
(SNR) and the constant £m set the accuracy level. We 
set in this work em = 0.005 or em = 0.035 as suggested 
in [l9[ (see references therein). 

Conditions (|A4[) and (|A5|) depend on the distance of 
the detector form the source. If they are written in term 
of the inaccuracy functional, the dependence on absolute 
scales can be moved to the right-hand-side and expressed 
only in term of the SNR. The accuracy requirements then 
read, 




faithful, 
effectual, 



(A6) 



where the functional dependence has been omitted for 
clarity. The level e < 1 is here set as e = 0.5. Equivalent 
conditions to Eq. (|A6[) can be expressed in term of the L2 
norm of the time domain signals (l9| . however they are 
not considered here because they seem more restrictive 
than the frequency domain criteria. 

In order to evaluate the accuracy of NR waveforms, 
one can consider h x as the best waveform model (the 
extrapolated in resolution in our case), and construct 
h m = h x + Sh from the error estimates (the last resolu- 
tion waveform in our case). The accuracy requirements 
in Eq. (|A6|) then quantify the accuracy of the NR wave- 
forms. In the analysis presented in the paper the integra- 
tion interval / e [0, oo] is approximated as / € [/o, /max], 
covering only the inspiral physical frequencies. 

The Wiener product is computed by scaling the wave- 
forms to physical units and to an effective distance, D e g, 
typically given in 100 Mpc. From the code output rh 22 , 
we: (i) recover rh + from the (2,2) multipole, using 
the expression for the spin weighted spherical harmon- 
ics, - 2 Y 2 ± 2 (6,(j)) = vV(64tt) exp(±i2<£)(l ± cos(9) 2 ; and 
(ii) scale r /i + to an effective distance. Assuming the ra- 
diation is emitted on the z axis, perpendicularly to the 
orbital plane, one has, 



r ^'- 2v u 1 - 2 



rh+(t) 

h + (t,D cff ) = rh+{t) GM & c 



Y 22 h 22 + - z Y 2 _ 2 h 2 _ 2 ) (A7) 
0.6308 r SR (h 22 ) (for 9 = , = 0) 
_ 2 /Dcff y 1 



\Mpc ) 
rh + (t) 4.7857 x 10~ 20 



(A8) 



/ Deff V 

\MpcJ 



Unit conversion: 1 Mpc a 
GMqc' 2 ~ 1.47670133 x 10 5 



3.08568025 x 10 24 cm, 
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